/*------------------------------------------------------------------------
 *  Copyright (C) 2010  Luis M. de la Cruz
 *
 *  This file is part of TUNA
 *
 *  TUNA is free software: you can redistribute it and/or modify
 *  it under the terms of the GNU General Public License as published by
 *  the Free Software Foundation, either version 3 of the License, or
 *  (at your option) any later version.
 *
 *  TUNA is distributed in the hope that it will be useful,
 *  but WITHOUT ANY WARRANTY; without even the implied warranty of
 *  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
 *  GNU General Public License for more details.
 *
 *  You should have received a copy of the GNU General Public License
 *  along with this program.  If not, see <http://www.gnu.org/licenses/>.
 ------------------------------------------------------------------------*/

namespace Tuna {

  /*!
   * Gauss-Seidel solution from problems in 1D.
   */
  template<typename Teq>
  int Solver::GaussSeidel1D(Teq &eq, double tolerance, int max_iter)
  {
    typedef typename Typeinfo<Teq>::prec_t prec_t;  
    prec_t residual = tolerance + 1.0;
    int count_iter = 0;
    static const int bi = eq.get_bi();
    static const int ei = eq.get_ei();
    Range I(bi, ei);
    
    while( (residual > tolerance) && (count_iter < max_iter) ) {
      eq.phi(I) = ( eq.aE(I) * eq.phi(I+1) + eq.aW(I) * eq.phi(I-1) + eq.sp(I)) / eq.aP(I);
      residual = eq.calcResidual();
      //      std::cout << "\n" << count_iter << "\t Residual = " << residual;
      count_iter++;
    }
    return count_iter;
  }
  
/*!
 ***************************************************************************
 * This function arranges the arrays of the Tri-diagonal linear system that
 * comes from 1D equations, in the manner required by the Thomas algorithm.
 * Then it solves the system using that algorithm. 
 ***************************************************************************
 *  \author Luis M. de la Cruz [ Mon May  4 21:31:07 CDT 2009 ]
 ***************************************************************************
 */
  template<class Teq>
  bool Solver::TDMA1D(Teq &eq) 
  {
    static const int bi = eq.get_bi();
    static const int ei = eq.get_ei();
    Range I(bi, ei);
    eq.phi(I) = Solver::thomas( eq.aP(I),   // a - main diagonal
				eq.aE(I),   // b - upper diagonal
				eq.aW(I),   // c - lower diagonal
				eq.sp(I) ); // d - left-side
    return 0;
  }

  /*!
   * In this function solves the system of equations until a tolerance is
   * achieved or a number of iteration is reached. The functions lineX_2D
   * and lineY_2D are called alternately to obtain better results. This
   * function works for 2D problems.
   */
  template<class Teq>
  int Solver::TDMA2DX(Teq &eq, double tolerance, int max_iter, double theta = 1.85) 
  {   
    typedef typename Typeinfo<Teq>::prec_t prec_t;  
    prec_t residual = tolerance + 1.0;
    int count_iter = 0;
  
    while( (residual > tolerance) && (count_iter < max_iter) ) {
      lineX_2D(eq, theta); 
      lineY_2D(eq, theta); 
      residual = eq.calcResidual();
      //      std::cout << "\n" << count_iter << "\t Residual = " << residual;
      count_iter++;
    }
    return count_iter;
  }

  /*!
   * In this function solves the system of equations until a tolerance is
   * achieved or a number of iteration is reached. The functions lineY_2D
   * and lineX_2D are called alternately to obtain better results.  This
   * function works for 2D problems.
   */
  template<class Teq>
  int Solver::TDMA2DY(Teq &eq, double tolerance, int max_iter, double theta = 1.85) 
  {   
    typedef typename Typeinfo<Teq>::prec_t prec_t;  
    prec_t residual = tolerance + 1.0;
    int count_iter = 0;
        
    while( (residual > tolerance) && (count_iter < max_iter) ) {
      lineY_2D(eq, theta); 
      lineX_2D(eq, theta); 
      residual = eq.calcResidual();
      //      std::cout << "\n" << count_iter << "\t Residual = " << residual;
      count_iter++;
    }
    return count_iter;
  }

  /*!
   * In this function solves the system of equations until a tolerance is
   * achieved or a number of iteration is reached. The functions lineY_3D,
   * lineZ_3D and lineX_3D are called alternately to obtain better results. 
   * This function works for 3D problems.
   */
  template<class Teq>
  int Solver::TDMA3D(Teq &eq, double tolerance, int max_iter, double theta = 1.85) 
  {
    typedef typename Typeinfo<Teq>::prec_t prec_t;  
    prec_t residual = tolerance + 1.0;
    int count_iter = 0;

    while( (residual > tolerance) && (count_iter < max_iter) ) {
      lineY_3D(eq, theta);
      lineZ_3D(eq, theta);
      lineX_3D(eq, theta);
      residual = eq.calcResidual();
      //      std::cout << "\n" << count_iter << "\t Residual = " << residual;
      count_iter++;
    }
    return count_iter;
  }

  /*!
   * This function arranges the arrays of a 2D mesh, in order to be passed to
   * the Thomas algorithm. Once the arrays are formed, the Thomas algorithm
   * is applied to find the solution in one line of the mesh in the x-direction.
   * Then we sweept all the lines moving on the y-direction.
   */
template<class Teq>
bool Solver::lineX_2D(Teq &eq, double theta) {
//
// VAN DOORMAL AND RAITHBY (84) : theta = 1 --> original TDMA
//
    typedef typename Typeinfo<Teq>::prec_t prec_t; 
    typedef typename TunaArray<prec_t, 1>::huge array;

    static const int bi = eq.get_bi();
    static const int ei = eq.get_ei();
    static const int bj = eq.get_bj();
    static const int ej = eq.get_ej();
    Range I(bi, ei); 

    int N = ei - bi + 1;
    static array a(N), d(N);

    for(int j = bj; j <= ej; ++j) {
      for(int i = bi, count = 0; i <= ei; ++i, ++count) {
	a(count) = eq.aP(i,j) - eq.aN(i,j) * (theta - 1);
	d(count) = eq.sp(i,j) + 
	  eq.aN(i,j) * (eq.phi(i,j+1) - eq.phi(i,j) * (theta - 1)) + 
	  eq.aS(i,j) * eq.phi(i,j-1);
      }
      eq.phi(I, j) = Solver::thomas( a, eq.aE(I,j), eq.aW(I,j) ,d );
    }
    return 0;
}

  /*!
   * This function arranges the arrays of a 2D mesh, in order to be passed to
   * the Thomas algorithm. Once the arrays are formed, the Thomas algorithm
   * is applied to find the solution in one line of the mesh in the y-direction.
   * Then we sweept all the lines moving on the x-direction.
   */
template<class Teq>
bool Solver::lineY_2D(Teq &eq, double theta) {
//
// VAN DOORMAL AND RAITHBY (84) : theta = 1 --> original TDMA
//
   typedef typename Typeinfo<Teq>::prec_t prec_t; 
   typedef typename TunaArray<prec_t, 1>::huge array;

    static const int bi = eq.get_bi();
    static const int ei = eq.get_ei();
    static const int bj = eq.get_bj();
    static const int ej = eq.get_ej();
    Range J(bj, ej); 

    int N = ej - bj + 1;
    static array a(N), d(N);

    int pausa;

    for(int i = bi; i <= ei; ++i) {
	for(int j = bj, count = 0; j <= ej; ++j, ++count) {
	    a(count) = eq.aP(i,j) - eq.aE(i,j) * (theta - 1);
	    d(count) = eq.sp(i,j) + 
		eq.aE(i,j) * (eq.phi(i+1,j) - eq.phi(i,j) * (theta - 1)) +
		eq.aW(i,j) * eq.phi(i-1,j);
	}

	if (BANDERA == 1 && i == bi) {
	  cout << "\n PHI = " << a; 
	  cout << "\n PHI = " << eq.aN(i,J); 
	  cout << "\n PHI = " << eq.aS(i,J); 
	  cout << "\n PHI = " << a; 
	  cin >> pausa;
	}
	eq.phi(i,J) = Solver::thomas( a, eq.aN(i,J), eq.aS(i,J), d );
    }

    return 0;
}

  /*!
   * This function arranges the arrays of a 3D mesh, in order to be passed to
   * the Thomas algorithm. Once the arrays are formed, the Thomas algorithm
   * is applied to find the solution in one line of the mesh in the x-direction.
   * Then we sweept all the lines in plane X-Z moving on the z-direction. Once
   * a plane is completed, we move upwards in y-direction and then we repeat the
   * process until the total of planes are covered.
   */
template<typename Teq>
bool Solver::lineX_3D(Teq &eq, double theta) {
//
// VAN DOORMAL AND RAITHBY (84) : theta = 1 --> original TDMA
//
    typedef typename Typeinfo<Teq>::prec_t prec_t; 
    typedef typename TunaArray<prec_t, 1>::huge array;

    static const int bi = eq.get_bi();
    static const int ei = eq.get_ei();
    static const int bj = eq.get_bj();
    static const int ej = eq.get_ej();
    static const int bk = eq.get_bk();
    static const int ek = eq.get_ek();
    Range I(bi, ei); 

    int N = ei - bi + 1;
    static array a(N), d(N);
    
    for(int j = bj; j <= ej; ++j) 
	for(int k = bk; k <= ek; ++k)
	{
	    for(int i = bi, count = 0; i <= ei; ++i, ++count) {
	      a(count) = eq.aP(i,j,k) - eq.aN(i,j,k) * (theta - 1);
	      d(count) = eq.sp(i,j,k) + 
		    eq.aN(i,j,k) * (eq.phi(i,j+1,k) -
				    eq.phi(i,j,k) * (theta - 1)) + 
		    eq.aS(i,j,k) * eq.phi(i,j-1,k) + 
		    eq.aF(i,j,k) * eq.phi(i,j,k+1) +
		    eq.aB(i,j,k) * eq.phi(i,j,k-1);
	    }
	    eq.phi(I,j,k) = Solver::thomas( a, eq.aE(I,j,k), eq.aW(I,j,k) ,d );
	}
    return 1;
}

  /*!
   * This function arranges the arrays of a 3D mesh, in order to be passed to
   * the Thomas algorithm. Once the arrays are formed, the Thomas algorithm
   * is applied to find the solution in one line of the mesh in the y-direction.
   * Then we sweept all the lines in plane Y-Z moving on the z-direction. Once
   * a plane is completed, we move in x-direction and then we repeat the
   * process until the total of planes are covered.
   */
template<typename Teq>
bool Solver::lineY_3D(Teq &eq, double theta) {
//
// VAN DOORMAL AND RAITHBY (84) : theta = 1 --> original TDMA
//
    typedef typename Typeinfo<Teq>::prec_t prec_t; 
    typedef typename TunaArray<prec_t, 1>::huge array;

    static const int bi = eq.get_bi();
    static const int ei = eq.get_ei();
    static const int bj = eq.get_bj();
    static const int ej = eq.get_ej();
    static const int bk = eq.get_bk();
    static const int ek = eq.get_ek();
    Range J(bj,ej);

    int N = ej - bj + 1;
    static array a(N), d(N);

    for(int i = bi; i <= ei; ++i) 
	for(int k = bk; k <= ek; ++k)
	{
	    for(int j = bj, count = 0; j <= ej; ++j, ++count) {
		a(count) = eq.aP(i,j,k) - eq.aE(i,j,k) * (theta - 1);
		d(count) = eq.sp(i,j,k) + 
		    eq.aE(i,j,k) * (eq.phi(i+1,j,k) - 
				    eq.phi(i,j,k) * (theta - 1)) +
		    eq.aW(i,j,k) * eq.phi(i-1,j,k) +
		    eq.aF(i,j,k) * eq.phi(i,j,k+1) +
		    eq.aB(i,j,k) * eq.phi(i,j,k-1);
	    }
	    eq.phi(i,J,k) = Solver::thomas( a, eq.aN(i,J,k), eq.aS(i,J,k), d );
	}
    
    return 1;
}

  /*!
   * This function arranges the arrays of a 3D mesh, in order to be passed to
   * the Thomas algorithm. Once the arrays are formed, the Thomas algorithm
   * is applied to find the solution in one line of the mesh in the z-direction.
   * Then we sweept all the lines in plane Y-Z moving on the y-direction. Once
   * a plane is completed, we move in x-direction and then we repeat the
   * process until the total of planes are covered.
   */
template<typename Teq>
bool Solver::lineZ_3D(Teq &eq, double theta) {
//
// VAN DOORMAL AND RAITHBY (84) : theta = 1 --> original TDMA
//
    typedef typename Typeinfo<Teq>::prec_t prec_t; 
    typedef typename TunaArray<prec_t, 1>::huge array;

    static const int bi = eq.get_bi();
    static const int ei = eq.get_ei();
    static const int bj = eq.get_bj();
    static const int ej = eq.get_ej();
    static const int bk = eq.get_bk();
    static const int ek = eq.get_ek();
    Range K(bk,ek);

    int N = ek - bk + 1;
    static array a(N), d(N);

    for(int i = bi; i <= ei; ++i) 
	for(int j = bj; j <= ej; ++j) 
	{
	    for(int k = bk, count = 0; k <= ek; ++k, ++count) {
		a(count) = eq.aP(i,j,k) - eq.aE(i,j,k) * (theta - 1);
		d(count) = eq.sp(i,j,k) + 
		    eq.aE(i,j,k) * (eq.phi(i+1,j,k) - 
				    eq.phi(i,j,k) * (theta - 1)) +
		    eq.aW(i,j,k) * eq.phi(i-1,j,k) +
		    eq.aN(i,j,k) * eq.phi(i,j+1,k) +
		    eq.aS(i,j,k) * eq.phi(i,j-1,k);
	    }
	    eq.phi(i,j,K) = Solver::thomas( a, eq.aF(i,j,K), eq.aB(i,j,K), d );
	}
    return 1;
}

/*!
 ***************************************************************************
 * Thomas algorithm for solving tridiagonal systems.
 * \verbatim
      |a1   b1    0    0    0  ...    0 |   | x1 |   | d1 |
      |c2   a2   b2    0    0  ...    0 |   | x2 |   | d2 |
      | 0   c3   a3   b3    0  ...    0 |   | x3 |   | d3 |
      | 0    0   c4   a4   b4  ...    0 | * | x4 | = | d4 |
      | .    .    .    .    .  ...    0 |   |  . |   |  . |
      | .    .    .    .    .  ...    0 |   |  . |   |  . |
      | 0    0    0    0  cN-1 aN-1 bN-1|   |xN-1|   |dN-1|
      | 0    0    0    0    0   cN   aN |   | xN |   | dN | \endverbatim
 *  
 * Step 1. \f$ P_1 = \frac{b_1}{a_1} \f$ and \f$ Q_1 = \frac{d_1}{a_1} \f$
 *
 * Step 2. For k = 2 to N do     (Forward step)
 * \f[
 *     P_k = \frac{b_k}{a_k - c_k P_{k-1}}
 *     Q_k = \frac{d_k - c_k Q_{k-1}}{a_k - c_k P_{k-1}} \f]
 *
 * Step 3. \f$ x_N = Q_N \f$
 *
 * Step 4. For k = N-1 to 1      (Backward substitution)
 * \f[ x_k = Q_k - P_k x_{k+1} \f]
 *
 * x contains the final solution.
 *
 * \param a diagonal of the matrix.
 * \param b upper diagonal
 * \param c lower diagonal
 * \param d left-side of the system
 ***************************************************************************
 *  \author  Luis M. de la Cruz [ Sat Mar 22 18:15:46 GMT 2008 ]
 ***************************************************************************/
template<typename Vector>
Vector Solver::thomas(const Vector &a, const Vector &b, 
		      const Vector &c, const Vector &d)
{
  int N = a.size();
  // The static doesn't work with the twophase problem the
  // the number of points is different on each axis ????
  //  static 
  Vector P(N), Q(N), x(N);
    
  P(0) = b(0) / a(0);  // Arrays begin in 0 in
  Q(0) = d(0) / a(0);  // this implementation

  double divisor;  
  for(int k = 1; k < N; ++k)     // Forward step
    {
      divisor = a(k) - c(k) * P(k-1);
      P(k) = b(k) / divisor;
      Q(k) = ( d(k) + c(k) * Q(k-1) ) / divisor;
    }
  x(N-1) = Q(N-1);    

  for(int k = N-2; k >= 0; --k)  // Backward substitution
    x(k) = P(k) * x(k+1) + Q(k);

  return x;
}

} // Tuna namespace



/*****  REVERSE SWAPPING
template<class T_eq>
bool Solver::lineX_2D_R(T_eq &eq) {
//
// VAN DOORMAL AND RAITHBY (84) : theta = 1 --> original TDMA
//
    typedef typename T_eq::array array;
    int bi = eq.bi, ei = eq.ei, bj = eq.bj, ej = eq.ej;
    int N = ei - bi + 1;
    array a(N), d(N);
    Range I(bi,ei);

    for(int j = ej; j >= bj; --j) {
	for(int i = bi, count = 0; i <= ei; ++i, ++count) {
	    a(count) = eq.aP(i,j) - eq.aS(i,j) * (eq.theta - 1);
	    d(count) = eq.sp(i,j) + 
		eq.aS(i,j) * (eq.phi(i,j-1) - eq.phi(i,j) * (eq.theta - 1)) + 
		eq.aN(i,j) * eq.phi(i,j+1);
	}
	eq.phi(I, j) = Solver::thomas( a, eq.aE(I,j), eq.aW(I,j) ,d );
    }
    
    return 1;
}

template<class T_eq>
bool Solver::lineY_2D_R(T_eq &eq) {
//
// VAN DOORMAL AND RAITHBY (84) : theta = 1 --> original TDMA
//
    typedef typename T_eq::array array;
    int bi = eq.bi, ei = eq.ei, bj = eq.bj, ej = eq.ej;
    int N = ej - bj + 1;
    array a(N), d(N);
    Range J(bj,ej);

    for(int i = ei; i >= bi; --i) {
	for(int j = bj, count = 0; j <= ej; ++j, ++count) {
	    a(count) = eq.aP(i,j) - eq.aW(i,j) * (eq.theta - 1);
	    d(count) = eq.sp(i,j) + 
		eq.aW(i,j) * (eq.phi(i-1,j) - eq.phi(i,j) * (eq.theta - 1)) +
		eq.aE(i,j) * eq.phi(i+1,j);
	}
	eq.phi(i,J).reverse(secondDim) = Solver::thomas(a.reverse(secondDim),
                                    eq.aN(i,J).reverse(secondDim),	       
                                    eq.aS(i,J).reverse(secondDim),
                                    d.reverse(secondDim) );
    }

    return 1;
}  
*****/








